
/*
 * En este archivo, implementaremos las funciones definidas
 * en el archivo matrices.h
 */

#include <stdlib.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include "matrices.h"


/*
 * Implementación de la función definida en cabecera.
 */
struct MATRIZ* crearMatriz(int dimension) {

    // No podemos crear matrices con dimension menor a 1.
    if (dimension<1) {
        return NULL;
    }

    // Creamos el vector de punteros a las filas
    double **matrizTmp;
    matrizTmp = (double **)malloc(sizeof(double *) * dimension);

    // Creamos una a una las filas
    int i;
    for (i=0;i<dimension;i++) {
        matrizTmp[i] = malloc(sizeof(double) * dimension);
    }

    // Creamos el objeto struct MATRIZ que retornaremos
    struct MATRIZ *matriz = malloc(sizeof(struct MATRIZ));
    matriz->dimension = dimension;
    matriz->elementos = matrizTmp;
    matriz->nombre = NULL;

    return matriz;
    
}


/*
 * Implementación de la función definida en cabecera.
 */
void eliminarMatriz(struct MATRIZ *matriz) {
    
    // Eliminamos la memoria ocupada por los elementos de la matriz
    int i;
    for (i=0; i < matriz->dimension ; i++) {
        double *fila = (matriz->elementos)[i];
        free(fila);
        fila = NULL;
    }
    
    // Eliminamos el array de punteros
    free(matriz->elementos);
    matriz->elementos = NULL;

    // Eliminamos el nombre
    free(matriz->nombre);
    matriz->nombre = NULL;

    // Eliminamos el objeto matriz
    free(matriz);
    matriz = NULL;

}


/*
 * Implementación de la función definida en cabecera.
 */
struct VECTOR* crearVector(int dimension) {
    
    // No existen vectores de dimension 0 o negativa!
    if (dimension<1)
        return NULL;
    
    // Reservamos la memoria necesaria para guardar el vector en memoria
    double *vectorTmp = (double *)malloc(sizeof(double)*dimension);

    // Creamos el objeto que devolveremos y lo devolvemos
    struct VECTOR * vector = malloc(sizeof(struct VECTOR) * dimension);
    vector->dimension = dimension;
    vector->elementos = vectorTmp;
    vector->nombre = NULL;
    
    return vector;
    
}


/*
 * Implementación de la función definida en cabecera.
 */
void eliminarVector(struct VECTOR *vector) {
    
    // Eliminamos el vector de elementos
    free(vector->elementos);
    vector->elementos = NULL;

    // Eliminamos el nombre
    free(vector->nombre);
    vector->nombre = NULL;

    // Eliminamos el objeto vector
    free(vector);
    vector = NULL;
    
}


/*
 * Implementación de la función definida en cabecera.
 */
void multiplicar(struct MATRIZ *matriz, struct VECTOR *vector, struct VECTOR *resultado) {
    
    // Declaracion de variables
    int i,j;            // Indices de iteracion
    double suma;        // Buffer para la suma en la multiplicacion de filas por columnas

    // Calculamos cada elemento del resultado
    for(i=0; i < resultado->dimension ; i++) {
        
        suma = 0; // Reiniciamos el buffer
        
        // Hacemos la sumatoria de la multiplicaciones de los elementos
        for(j=0; j < resultado->dimension ; j++) {
            suma += matriz->elementos[i][j] * vector->elementos[j];
        }

        // Asignamos la sumatoria al elemento correspondiente del resultado
        resultado->elementos[i] = suma;
        
    }

}


/*
 * Implementación de la función definida en cabecera.
 */
struct MATRIZ* multiplicarMatricesLU(struct MATRIZ *matrixU, struct MATRIZ *matrizL) {

    // Declaracion de variables
    struct MATRIZ *resultado = crearMatriz(matrizL->dimension);
    int tam = matrizL->dimension;
    int i,j;                // Indices de iteracion
    register int k;         // Indices de interacion optimizado
    double suma;   // Buffer para la suma en la multiplicacion de filas por columnas
    
    // Buffereamos las indirecciones para acelerar el cálculo
    double *U_i;
    double *resultado_i;
    
    int limite; // Aqui guardaremos hasta donde iteraremos el k (evitando multplicaciones por 0)

    // Buffereamos los punteros asi evitamos desreferencias.
    double **U_elem = matrixU->elementos;
    double **L_elem = matrizL->elementos;
    double **resultado_elem = resultado->elementos;
    
    for(i=0 ; i < tam ; i++) {
        
        // Buffereamos una indireccion en ambas variables, por motivos de rendimiento
        U_i = U_elem[i];
        resultado_i = resultado_elem[i];
        
        for(j=0 ; j < tam ; j++) {
            
            limite = (i<=j) ? i : j ;
            suma = 0;
            
            for(k=0;k<tam;k++) {
                suma += U_i[k] * L_elem[k][j];
            }
            
            resultado_i[j] = suma;
            
        }
    }

    // Devolvemos la matriz resultante
    return resultado;
        
}
